%mimo mmse V_BLAST bpsk rayleigh channel
%modulation:bspk
%channel coding:NULL
%mimo channel:rayleigh
%detection method : MMSE-SIC ZF-SIC ZF MMSE

clear all;
N = 100000; % number of transmit bits N 为Tx的倍数
snr_dB = [0:15]; % Eb/No values
Tx = 2;
Rx = 2;

err_all1 = [];
err_all2 = [];
err_all3 = [];
err_all4 = [];

for ii = 1:length(snr_dB)
    snr_val = sqrt(10^(-snr_dB(ii)/10));
    source = round(rand(1,N));
    ssmod = 2*source - 1; % bpsk modulation 0->-1 1->1
    ssmod = repmat(ssmod,Rx,1);
    ssmod = reshape(ssmod, Rx, Tx, N/Tx);

    % 每根天线上数据的长度为N/Tx

    %generate the channle
    h = 1/sqrt(2)*[randn(Rx,Tx, N/Tx) + j*randn(Rx, Tx, N/Tx)]; % rayleigh channel
    n = 1/sqrt(2)*[randn(Rx, N/Tx) + j*randn(Rx, N/Tx)];        % awgn 0dB variance

    % received signal
    y = reshape(sum(h.*ssmod,2),Rx, N/Tx, 1)+ snr_val*n; %size : Rx*(N/Tx)
%      y = reshape(sum(h.*ssmod,2),Rx, N/Tx, 1);


    % receiver detection
    % --------------
    % Forming the MMSE equalization matrix w = inv(H'*H + sigma^2)*H'
    % inverse of a [2*2] matrxi  [a b; c d] = 1/(ad-bc)[d -b;-c a]

    for jj = 1:N/Tx
        for t = 1:Tx

            % --------- ZF detection
            g1 = pinv(h(:,:,jj)); % 对信道矩阵求逆
            y_temp1(:,jj) = g1*y(:,jj);
            yy1(:,jj) = real(y_temp1(:,jj))>0;
            % -------- MMSE detection
            I = diag(ones(1,Tx));
            hj = h(:,:,jj);
            g2 = inv(hj'*hj + (snr_val^2)*I)*hj';
            y_temp2(:,jj) = g2*y(:,jj);
            yy2(:,jj) = real(y_temp2(:,jj)) > 0;

        end
    end

    for n = 1:N/Tx
        set1 = zeros(Tx,1);
        set2 = zeros(Tx,1);
        hn = h(:,:,n);
        yr_temp1 = y(:,n);
        yr_temp2 = y(:,n);
        hn_temp1 = hn;
        hn_temp2 = hn;
        for r = 1:Tx
            % ---------- ZF-SIC detection
          
            w1 = pinv(hn_temp1);
            fw1 = sum(abs(w1).^2,2)+set1; % 每一行的二范数
      

            [Tmin, k1] = min(fw1); % 找出最小的行

            ww1 = w1(k1,:);
            d_temp1 = ww1*yr_temp1;
            dd1(k1,n) = real(d_temp1) > 0; % QPSK hard decisio
    
            yr_temp1 = yr_temp1 -(2 * dd1(k1,n) - 1)*hn_temp1(:,k1);
            hn_temp1(:,k1) = zeros(Rx, 1); % 把H的第K列置为0
            set1(k1,1) =inf;

            % ---------- MMSE-SIC detection
           
            I = diag(ones(1,Tx));
            w2 = inv(hn_temp2'*hn_temp2 + (snr_val^2)*I)*hn_temp2';
            fw2 = sum(abs(w2).^2,2)+set2; % 每一行的二范数
%             fw2(set2) = inf;

            k2 = find(fw2 == min(fw2)); % 找出最小的行

            ww2 = w2(k2,:);
            d_temp2 = ww2*yr_temp2;
            dd2(k2,n) = real(d_temp2) > 0; % QPSK hard decision

            yr_temp2 = yr_temp2 - hn_temp2(:,k2)*(2 * dd2(k2,n) - 1);
            hn_temp2(:,k2) = zeros(Rx, 1); % 把H的第K列置为0
            set2(k2) = inf;

        end
    end


    yz = reshape(yy1,1,numel(yy1)); % the output of ZF detection
    ym = reshape(yy2,1,numel(yy2)); % the output of MMSE detection
    ysz = reshape(dd1,1,numel(dd1)); % the output of ZF-SIC detection
    ysm = reshape(dd2,1,numel(dd2)); % the output of MMSE-SIC detection



    error1 = length(find(yz~=source))/N;
    error2 = length(find(ym~=source))/N;
    error3 = length(find(ysz~=source))/N;
    error4 = length(find(ysm~=source))/N;

    err_all1(ii) = error1;
    err_all2(ii) = error2;
    err_all3(ii) = error3;
    err_all4(ii) = error4;

end

save result.mat err_all1 err_all2 err_all3 err_all4

EbN0Lin = 10.^(snr_dB/10);
theoryBer_nRx1 = 0.5.*(1-1*(1+1./EbN0Lin).^(-0.5));
% close all;
semilogy(snr_dB,theoryBer_nRx1,'bp-','LineWidth',2);
hold on
semilogy(snr_dB, err_all1,'->');
semilogy(snr_dB, err_all2,'r-o');
semilogy(snr_dB, err_all3,'-*');
semilogy(snr_dB, err_all4,'r-d');


legend('theory(Tx=2,Rx=2, ZF)', 'sim(Tx=2, Rx=2, ZF', 'sim(Tx=2, Rx=2, MMSE)',...
    'sim (Tx=2, Rx=2, ZF-SIC)', 'sim(Tx=2, Rx=2, MMSE-SIC)')
xlabel('snr(dB)');
ylabel('BER');


% the end

